#1.TWAS.meta.p.value
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201204")
sel=list.files(pattern = "TWAS.BD")
library(meta)
f1=read.csv(sel[1],header=T)
f2=read.csv(sel[2],header=T)
f3=read.csv(sel[3],header=T)
info.data=data.frame(Term=c("117K","61K","53425","8544","807","200","3969"),
                     unique.gene.num=c(16331,13621,12613,4231,879,236,2308))

result=data.frame(matrix(,nrow = 7,ncol = 5))
names(result)=c("Term","OR","P.value","upper","lower")
for(i in 1:7){
case.overlap.num=as.numeric(rep(f1[i+1,"overlap.num"],3))
case.total.num=as.numeric(rep(info.data[i,"unique.gene.num"],3))
con.overlap.num=as.numeric(c(f1[1,"overlap.num"],f2[1,"overlap.num"],f3[1,"overlap.num"]))
con.total.num=as.numeric(c(f1[1,"OR"],f2[1,"OR"],f3[1,"OR"]))

tdata=data.frame(case.overlap.num,case.total.num,con.overlap.num,con.total.num)
metaor3<-metabin(case.overlap.num,case.total.num,con.overlap.num,con.total.num,data=tdata,sm="OR")
result[i,1]=info.data[i,1]
result[i,2]=OR=exp(metaor3$TE.fixed)
result[i,3]=metaor3$pval.fixed
result[i,4]=upper=exp(metaor3$upper.fixed)
result[i,5]=lower=exp(metaor3$lower.fixed)
}
result$FDR=p.adjust(result$P.value,method = "BH")
write.csv(result,"./meta.p.TWAS.BD.statis.csv",quote=F,row.names = F)

#2.AJHG
library(meta)
sel=list.files(pattern = "AJHG")
f1=read.csv(sel[1],header=T)
f2=read.csv(sel[2],header=T)
f3=read.csv(sel[3],header=T)
info.data=data.frame(Term=c("117K","61K","53425","8544"),
                     unique.num=c(117012,61725,53425,8544))

result=data.frame(matrix(,nrow = 4,ncol = 5))
names(result)=c("Term","OR","P.value","upper","lower")
for(i in 1:4){
case.overlap.num=as.numeric(rep(f1[i+1,"overlap.num"],3))
case.total.num=as.numeric(rep(info.data[i,"unique.num"],3))
con.overlap.num=as.numeric(c(f1[1,"overlap.num"],f2[1,"overlap.num"],f3[1,"overlap.num"]))
con.total.num=as.numeric(rep(150000,3))

tdata=data.frame(case.overlap.num,case.total.num,con.overlap.num,con.total.num)
metaor3<-metabin(case.overlap.num,case.total.num,con.overlap.num,con.total.num,data=tdata,sm="OR")
result[i,1]=info.data[i,1]
result[i,2]=OR=exp(metaor3$TE.fixed)
result[i,3]=metaor3$pval.fixed
result[i,4]=upper=exp(metaor3$upper.fixed)
result[i,5]=lower=exp(metaor3$lower.fixed)
}
result$FDR=p.adjust(result$P.value,method = "BH")
write.csv(result,"./meta.p.AJHG.statis.csv",quote=F,row.names = F)

#3.220K ASM
library(meta)
sel=list.files(pattern = "220K")
f1=read.csv(sel[1],header=T)
f2=read.csv(sel[2],header=T)
f3=read.csv(sel[3],header=T)
info.data=data.frame(Term=c("117K","61K","53425","8544","807","200"),
                     unique.num=c(117012,61725,53425,8544,807,200))

result=data.frame(matrix(,nrow = 6,ncol = 5))
names(result)=c("Term","OR","P.value","upper","lower")
for(i in 1:6){
case.overlap.num=as.numeric(rep(f1[i+1,"overlap.num"],3))
case.total.num=as.numeric(rep(info.data[i,"unique.num"],3))
con.overlap.num=as.numeric(c(f1[1,"overlap.num"],f2[1,"overlap.num"],f3[1,"overlap.num"]))
con.total.num=as.numeric(rep(150000,3))

tdata=data.frame(case.overlap.num,case.total.num,con.overlap.num,con.total.num)
metaor3<-metabin(case.overlap.num,case.total.num,con.overlap.num,con.total.num,data=tdata,sm="OR")
result[i,1]=info.data[i,1]
result[i,2]=OR=exp(metaor3$TE.fixed)
result[i,3]=metaor3$pval.fixed
result[i,4]=upper=exp(metaor3$upper.fixed)
result[i,5]=lower=exp(metaor3$lower.fixed)
}
result$FDR=p.adjust(result$P.value,method = "BH")
write.csv(result,"./meta.p.220K.ASM.statis.csv",quote=F,row.names = F)

#4.imprint.locis
library(meta)
sel=list.files(pattern = "imprint")[2:4]
f1=read.csv(sel[1],header=T)
f2=read.csv(sel[2],header=T)
f3=read.csv(sel[3],header=T)
info.data=data.frame(Term=c("117K","61K","53425","8544","807","3969","200","1195"),
                     unique.num=c(117012,61725,53425,8544,807,3969,200,1195))

result=data.frame(matrix(,nrow = 8,ncol = 5))
names(result)=c("Term","OR","P.value","upper","lower")
for(i in 1:8){
case.overlap.num=as.numeric(rep(f1[i+1,"overlap.num"],3))
case.total.num=as.numeric(rep(info.data[i,"unique.num"],3))
con.overlap.num=as.numeric(c(f1[1,"overlap.num"],f2[1,"overlap.num"],f3[1,"overlap.num"]))
con.total.num=as.numeric(rep(150000,3))

tdata=data.frame(case.overlap.num,case.total.num,con.overlap.num,con.total.num)
metaor3<-metabin(case.overlap.num,case.total.num,con.overlap.num,con.total.num,data=tdata,sm="OR")
result[i,1]=info.data[i,1]
result[i,2]=OR=exp(metaor3$TE.fixed)
result[i,3]=metaor3$pval.fixed
result[i,4]=upper=exp(metaor3$upper.fixed)
result[i,5]=lower=exp(metaor3$lower.fixed)
}
result$FDR=p.adjust(result$P.value,method = "BH")
write.csv(result,"./meta.p.imprinted.locis.statis.csv",quote=F,row.names = F)